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Abstract. We describe the construction of a conserved reaction-diffusion system that exhibits self-organized 
critical (avalanche-like) behavior under the action of a slow addition of particles. The model provides an 
illustration of the general mechanism to generate self-organized criticality in conserving systems. Extensive 
simulations in d — 2 and 3 reveal critical exponents compatible with the universality class of the stochastic 
Manna sandpile model. 
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Since the introduction of the Bak, Tang, and Wiesen- 
feld sandpile model Q , the concept of self-organized criti- 
■ cality (SOC) Q has witnessed a real explosion of activity, 
' covering both the description of new models and the pro- 
posal of several theoretical approaches, aiming at an un- 
derstanding of SOC phenomena in terms of standard sta- 
tistical mechanical concepts. At this respect, it has been 
shown that SOC in sandpile models is related to the be- 
| havior of absorbing-state phase transitions (APT) with 
. many absorbing states Indeed, this very idea is un- 

derlying a recipe proposed for the construction of SOC 
models ^ : any cellular automata, defined with conserved 
microscopic rules, and possessing many absorbing states, 
will display SOC behavior if slowly driven with the addi- 
tion of energy/particles at a rate h and dissipation at a 
rate e: i.e., in the double limit h — > 0, e — > 0, with h/e — > 
H . This mechanism is easily seen at work in all standard 
sandpile models proposed so far |l],f?],|]]. 

Given that most SOC systems are defined in terms of 
sandpile-like models (with the exception of the forest-fire 
H and extremal JhJ models), it becomes all the most in- 
terestingto explore the possibility of applying the recipe 
' of Ref. H to models of a qualitatively different sort. In 
this paper we consider a reaction-diffusion (RD) model 
showing an APT that conserves the total number of parti- 
cles []Tl[p~2| . This model exhibits a non-equilibrium phase 
transition in the same universality class of fixed energy 
stochastic sandpiles |^,^2). Here, we show that implement- 
ing the slow driving condition, the model reaches a station- 
ary state with an avalanche-like reaction activity with crit- 
ical properties. By measuring usual magnitudes character- 
izing the SOC behavior, we compare the model with stan- 
dard slowly driven sandpiles. The critical exponents mea- 
sured confirm the shared universality class with stochastic 
sandpiles, and provide a vivid illustration of the SOC gen- 
erating mechanism 



We focus on the two species RD system jy],[L^, re- 
cently proposed to describe APT coupled to a non-diffusive 
conserved field |l3| . The RD system is defined by the fol- 
lowing set of reaction steps: 

B^A with rate ki, (1) 
B + A -> 2B with rate k 2 . (2) 

In this system, B particles diffuse with diffusion rate Db, 
and A particles do not diffuse, that is, = 0. From the 
rate Eqs. (||) and it is clear that the dynamics con- 
serves the total number of particles N — Na + Nb, where 
Ni is the number of particles i = A, B. In this model, the 
dynamics is exclusively due to B particles, that we iden- 
tify as active particles. A particles do not diffuse and can- 
not generate spontaneously B particles. More specifically, 
A particles can only move via the motion of B particles 
that later on transform into A through Eq. (||). This im- 
plies that any configuration devoid of B particles is an 
absorbing state in which the system is trapped forever. 
The number of these absorbing states is infinite — in the 
thermodynamic limit — corresponding to all the possible 
redistributions of N particles of type A in the system. 
This RD process exhibits a phase transition from an ac- 
tive phase (with an everlasting activity of B particles) to 
an absorbing phase (no B particles) for a critical value 
p = p c of the total particle density Q • 

Here, we define a driven-dissipative version of the RD 
model by applying the recipe of Refs. On hypercubic 

lattices of size L with open boundary conditions, each site 
i stores a number et^ of A particles and bi of B particles. 
The occupation numbers on and bi can have any integer 
value, including = bi = 0, that is, void sites with no 
particles. The model is thus representing the dynamics of 
bosonic particles. The initial configuration is constructed 
by randomly distributing a number Nq of A particles in 
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the lattice. The initial occupation numbers a,i have a pois- 
sonian distribution, while bi = 0,Vi. Any configuration is 
stable whenever it fulfills this condition, i.e., in the absence 
of B particles. The system is driven by adding one B par- 
ticle to a randomly chosen site. A state with at least one 
B particle is called active. Active states evolve in time ac- 
cording to the following update rules, that mimic the dif- 
fusion and reaction steps in the RD system: I) Diffusion: 
on each lattice site, each B particle moves into a randomly 
chosen nearest neighbor site with probability 2d/(2d+ 1), 
and stays in the same site with probability l/(2d + 1); 
this results in an effective diffusivity Db = l/(2d + 1). 
II) After all sites have been updated for diffusion, we per- 
form the reactions: a) On each lattice site, each B particle 
is turned into an A particle with probability r%. b) At 
the same time, each A particle becomes a B particle with 
probability 1 — (1 — r2) bi , where bi is the total number 
of B particles in that site. This corresponds to the aver- 
age probability for an A particle of being involved in the 
reaction (|^) with any of the B particles present on the 
same site. The probabilities r\ and T2 are related to the 
reaction rates fci and hi defined in Eqs. (|l|) and (||). In 
general, we have that ri(ki = 0) = 0, rj(fcj = oo) = 1, 
and Ti is an increasing function of fcj. The analytic ex- 
pression of Ti as a function of fcj is presumably quite com- 
plex and nontrivial. However, as we will argue later, the 
knowledge of the precise relationship between Ti and hi 
is not necessary, since the critical behavior of the model 
should be independent of the exact values of the parame- 
ters n selected. B particles on boundary sites may choose 
to diffuse out of the lattice. In this case, the particle is re- 
moved out of the system, contributing to the dissipation. 
The system is updated in parallel until there are no more 
B particles and it is again in an absorbing state. During 
the dynamic evolution, the addition of new B particles is 
suspended; this of course corresponds to the slow-driving 
condition. The sequence of updates in the system (from 
the time we introduce a new B particle until a stable state 
is reached) is interpreted as an avalanche. We characterize 
avalanches by their size s and their duration t. The size 
of an avalanche is defined as the number of B particles 
present in the system at each time step, summed over all 
the parallel updates required to reach a new stable state. 
The duration of an avalanche is defined as the total num- 
ber of parallel updates performed during the avalanche. 
In the slow driving perspective, the existence of a critical 
stationary state is easily understood. Particles are added 
only in the absence of activity (p < p c ), while dissipation 
acts only during activity (p > p c ). This implies that dtp 
always drives the system toward p c , that in the thermo- 
dynamic limit is the only possible stationary value of the 
density gg. 

We have performed numerical simulations of this model 
in dimensions 2 and 3, with system sizes ranging from 
L = 64 to L = 1024 in d = 2, and from L = 74 to 
L = 280 in d = 3. The reaction rates r, reported here 
are r\ — 0.3 and r-i = 0.4 in d — 2, and r% = 0.4 and 
T2 = 0.5 in d = 3; the numerical results however, are in- 
dependent of the particular values of the reaction rates 
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Fig. 1. Number of A particles Na (dashed line) and average 
avalanche size (s) (full line) as a function of the number of 
avalanches T in a slowly driven conserved reaction diffusion 
system in a lattice of size L — 256. The initial state is an 
empty lattice. 

Ti employed. The independence of the dynamics from the 
particular parameters rj chosen can be easily grasped by 
the mean- field approach reported in Ref . Ill] , that shows 
the existence of a single critical point. More precisely, this 
can be understood in terms of renormalization-group ar- 
guments: The critical behavior of the model is ruled by 
the fixed point toward which the parameters flow under 
an appropriate renormalization-group transformation [ pT[ . 
Thus, the critical exponents are independent of the initial 
values of the parameters, and depend only on the value 
of the unique fixed point. It is then natural that simula- 
tions performed with different parameters will yield the 
same steady state behavior. This fact is confirmed by the 
numerical simulations, which show differences only in the 
transient regime. After the transient initial regime (whose 
length depends on L, fj, and the initial concentration Nq 
of A particles), the system reaches a steady state in which 
the stable configurations have a constant average number 
Na of A particles and avalanches have a constant average 
size (s). As an example, in Fig. |l|, we plot the number 
of A particles and the average avalanche size as a func- 
tion of the number of avalanches T in a simulation with 
system size L = 256 in d — 2, starting from an empty lat- 
tice. In analogy with sandpiles and SOC phenomena, in 
the steady state, we compute the probability distribution 
of sizes P(s) and times P(t) for the reaction events. By 
assuming the usual finite-size scaling form (FSS) Q 

P( S ,L) = S -^(^), (3) 

p{t,L) = t- r *g(J^j, (4) 

we can define the standard critical exponents, r s , r t , D 
(the fractal dimension), and z (the dynamic critical ex- 
ponent), which characterize the universality class of the 
model. Averages are performed over at least 5 x 10 6 ava- 
lanches in d = 2. The distributions on d = 3 are consid- 
erably noisier; averages have thus been done in this case 
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over 10 7 avalanches. In order to check the plausibility of 
the FSS hypothesis (||) and (0), and compute the corre- 
sponding critical exponents, we have applied the moment 
analysis technique, introduced in Refs. fll5|j . We define the 
q-th moment of the avalanche size distribution on a lattice 
of size L as (s q ) L = Jdss q P(s, L). If the FSS hypothesis 
(|J) is valid in the asymptotic limit of large s, then the q-th 
moment has the following dependence on system size: 



»(«) 



(5) 



The exponent o~ s (g) = D(q + 1 — r s ) is computed as the 
slope of the log-log plot of (s q ) L as a function of L. For 
large enough values of q (i.e., away from the region where 
the integral in (JsJ) is dominated by its lower cut-off), one 
can compute the fractal dimension D as the slope of a s (q) 
as a function of q: D = da s {q) / dq. Once D is known we 
can estimate t s using the relation er s (l) = D(2 — t s ). The 
exponent er s (l), giving the scaling of the average avalanche 
size, can be estimated using a standard argument in sand- 
piles: a new injected particle B performs a diffusing Brow- 
nian motion and has to travel, on average, a distance of 
order L 2 before reaching the boundary. In the station- 
ary state, to each B particle drop must correspond, on 
average, a B particle diffusing out of the system. This 
implies that the average avalanche size corresponds to 
the average number of diffusion steps needed for a B 
particle to reach the boundary; i.e., (s) ~ L 2 , and thus 
cr s (l) = 2. Along the same lines we can obtain the mo- 
ments of the avalanche time distribution. In this case, 
{t q )L ~ L at ^ q \ with da t {q)/dq = z. Analogous consid- 
erations for small q apply also for the time moment anal- 
ysis. Then, the r t exponent can be found using the scaling 
relation z{2 — Tt) = at(l). 

We have computed the exponents r s , r f , D, and z from 
our data, using the moment analysis method. Our results 
are reported in Tables [l] and 0. The validity of these ex- 
ponents can be checked a posteriori by means of a data 
collapse analysis: If the FSS hypothesis of Eqs. (||) and 
is correct, then the plots of the distributions, under the 
rescaling s — > s/L D and P(s) — > P{s)L Dt " , and corre- 
spondingly t -> t/L z and P(t) -»■ P(t)L ZTt , should col- 
lapse onto the same universal function, for different values 



D 



Tt 



Conserved RD 1.28(1) 2.75(1) 1.50(2) 1.54(1) 
Manna 1.28(1) 2.76(1) 1.48(2) 1.55(1) 

Table 1. Critical exponents for the conserved RD and the 
stochastic Manna models in d — 2. Figures in parenthesis in- 
dicate the statistical uncertainty in the last digit. Manna ex- 
ponents from Refs. ||,|l|,|l|,|§. 



D 



Tt 



Conserved RD 
Manna 



1.42(1) 
1.41(1) 



3.36(1) 
3.36(1) 



1.80(2) 
1.78(2) 



1.77(1) 
1.76(1) 



Table 2. Critical exponents for the conserved RD and the 
stochastic Manna models in d = 3. Figures in parenthesis in- 
dicate the statistical uncertainty in the last digit. 
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Fig. 2. Data collapse analysis of the integrated avalanche size 
distribution for the conserved RD model in d = 2. System sizes 
L = 64, 128, 256, 512, and 1024. 
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Fig. 3. Data collapse analysis of the integrated avalanche time 
distribution for the conserved RD model in d = 2. System sizes 
L = 64, 128, 256, 512, and 1024. 



of L. Figs. Hand|| show the data collapse for the distribu- 
tions of sizes and times in d — 2, respectively. Analogous 
plots for the case d = 3 are presented in Figs. ^ and M. 
From these plots we conclude that our model fulfills the 
FSS hypothesis, and is in this sense critical, being charac- 
terized by the exponents reported in Tables [l] and 0. 

Once we have shown that the slowly driven conserved 
RD model exhibits SOC behavior, it is interesting to check 
whether it shares the same universality class with any 
other SOC models. Given its stochastic rules, the natural 
candidate for comparison is the stochastic Manna sandpile 
model [||,^6|. In Table |l| we quote the exponents for the 
Manna model in d = 2, whose value has been relatively 
well established by different sets of independent simula- 
tions 0,^|,^9[|2^j. The case rf = 3 has not been stud- 
ied so thoroughly, with the exception of the simulations of 
Liibcck |lq| . Thus, in order to compare our results with the 
conserved RD model, we have carried out large-scale sim- 
ulations of the rf = 3 stochastic Manna model. The expo- 
nents obtained are given in Table ||, while Figures || and 
plot the data collapse for sizes and times, respectively. Av- 
erages are performed over 10 7 non-zero avalanches, and 
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Fig. 4. Data collapse analysis of the integrated avalanche size 
distribution for the conserved RD model in d = 3. System sizes 
L = 74, 100, 140, 200, and 280. 




Fig. 5. Data collapse analysis of the integrated avalanche time 
distribution for the conserved RD model in d = 3. System sizes 
L = 100, 140, 200, and 280. 




Fig. 6. Data collapse analysis of the integrated avalanche size 
distribution for the stochastic Manna model in d = 3. System 
sizes L = 100, 140, 200, 280, and 400. 
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Fig. 7. Data collapse analysis of the integrated avalanche time 
distribution for the stochastic Manna model in d = 3. System 
sizes L = 100, 140, 200, 280, and 400. 



the system sizes considered range from L = 74 to L — 400 
(larger than those achieved by Liibeck Jig]). 

From the examination of Tables |l| and ^, we conclude 
that the present conserved RD model exhibits exponents 
fully compatible with those of the stochastic Manna sand- 
pile model in both d = 2 and d = 3. This coincidence of ex- 
ponents proves that both models belong to the same uni- 
versality class. This fact is altogether not surprising, since 
both models exhibit the same basic symmetries: isotropic 
diffusion of particles, stochastic conserved microscopic rules, 
and presence of infinitely many — in the thermodynamic 
limit — absorbing states. This result represents a further 
confirmation of the universality claim made in Refs. |i2| , 
|l3| for this kind of systems. 

In summary, we have shown how to construct a slowly- 
driven conserved RD system which exhibits SOC behav- 
ior (avalanche macroscopic dynamics). The key points are 
the application of the slow-driving condition (addition of 
new B particles) plus boundary dissipation (diffusion of 
B particles out of the system). The model is character- 
ized by critical exponents that place it in the same uni- 
versality class than the Manna stochastic sandpile model. 
The related model proposed in Ref. |11 with Da ^ 0, 



however, belongs to a different universality class j^). The 
limit Da —* in the theory with Da ^ is non-analytic; 
any infinitesimal amount of diffusion in the energy field 
renormalizes to a finite value, and definitely changes the 
universality class of the model. 

This work has been supported by the European Net- 
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acknowledges support from the grant CICYT PB97-0693. 



References 

1. P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59, 
381 (1987). 

2. H. J. Jensen, Self- Organized Criticality (Cambridge Uni- 
versity Press, Cambridge, 1998). 

3. R. Dickman, A. Vespignani, and S. Zapperi, Phys. Rev. E 
57, 5095 (1998); A. Vespignani, R. Dickman, M. A. Mufioz, 
and S. Zapperi, Phys. Rev. Lett. 81, 5676 (1998). 

4. A. Vespignani, R. Dickman, M. A. Munoz, and S. Zapperi, 
Phys. Rev. E 62, 4564 (2000). 

5. R. Dickman, M. A. Munoz, A. Vespignani, and S. Zapperi, 
Braz. J. Phys. 30, 27 (2000). 



Romualdo Pastor-Satorras, Alessandro Vespignani: Reaction-diffusion system with self-organized critical behavior 



5 



6. G. Grinstein, in Scale Invariance, Interfaces and Nonequi- 
librium Dynamics, NATO Advanced Study Institute, Series 
B: Physics, vol. 344, A. McKane et al., Eds. (Plenum, New 
York, 1995). 

7. Y.-C. Zhang, Phys. Rev. Lett. 63, 470 (1989). 

8. S. S. Manna, J. Phys. A 24, L363 (1991). 

9. B. Drossel and F. Schwabl, Phys. Rev. Lett. 69, 1629 
(1992). 

10. P. Bak and K. Sneppen, Phys. Rev. Lett. 71, 4083 (1993). 

11. F. van Wijland, K. Oerding, and H. J. Hilhorst, Physica 
A 251, 179 (1998). 

12. R. Pastor-Satorras and A. Vespignani, Phys. Rev. E 62, 
R5875 (2000). 

13. M. Rossi, R. Pastor-Satorras, and A. Vespignani, Phys. 
Rev. Lett. 85, 1803 (2000). 

14. Finite Size Scaling, Vol. 2 of Current Physics-Sources and 
Comments, edited by J. L. Cardy (North Holland, Ams- 
terdam, 1988). 

15. M. De Menech, A. L. Stella, and C. Tebaldi, Phys. Rev. E 
58, R2677 (1998); C. Tebaldi, M. De Menech, and A. L. 
Stella, Phys. Rev. Lett 83, 3952 (1999). 

16. D. Dhar, Physica A 263, 4 (1999). 

17. A. Chessa, A. Vespignani, and S. Zapperi, Comput. Phys. 
Commun. 121-122, 299 (1999). 

18. S. Lubeck, Phys. Rev. E 61, 204 (2000). 

19. K. Nakanishi and K. Sneppen, Phys. Rev. E 55, 4012 
(1997). 

20. E. Milshtein, O. Biham, and S. Solomon, Phys. Rev. E 58, 
303 (1998). 



